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Elastic wave propagation is studied in a heterogeneous 2-D medium consisting of an elastic 
matrix containing randomly distributed circular elastic inclusions. The aim of this study is 
to determine the effective wavenumbers when the incident wavelength is similar to the radius 
of the inclusions. A purely numerical methodology is presented, with which the limitations 
usually associated with low scatterer concentrations can be avoided. The elastodynamic equa- 
tions are integrated by a fourth-order time-domain numerical scheme. An immersed interface 
method is used to accurately discretize the interfaces on a Cartesian grid. The effective field 
is extracted from the simulated data, and signal-processing tools are used to obtain the com- 
plex effective wavenumbers. The numerical reference solution thus-obtained can be used to 
check the validity of multiple scattering analytical models. The method is applied to the case 
of concrete. A parametric study is performed on longitudinal and transverse incident plane 
waves at various scatterers concentrations. The phase velocities and attenuations determined 
numerically are compared with predictions obtained with multiple scattering models, such as 
the Independent Scattering Approximation model, the Waterman- Truell model, and the more 
recent Linton-Martin-Conoir-Norris model. 



Keywords: ultrasounds; multiple scattering; effective medium; homogenization; numerical 
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1. Introduction 

We consider the propagation of elastic waves across a medium containing randomly 
distributed circular inclusions, the size of which is similar to that of the wavelength. 
The effective field, which is obtained by averaging the fields in all the possible 
disordered configurations, corresponds to that of waves propagating in an effective 
homogeneous medium. 

There exists three possible approaches for obtaining the effective wavenumbers 
(and equivalently, the effective phase velocity and attenuation): 

• theoretical approach, based on multiple-scattering models such as the classical 
Waterman- Truell [l| and Foldy [2] models. It provides closed- form expressions 
useful in practical applications. The main assumption is that the scatterer con- 
centration is low, i.e. typically less than 10 % [3]. At higher concentrations, more 
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sophisticated models developed in acoustics by Linton and Martin 0, [H], and 
extended to elasto dynamics by Conoir and Norris are required. But, to 

our knowledge, no rigorous error estimate is available, and the limits of validity 
are not accurately known; 
• experimental approach It introduces no limitation about the concentration 
of scatterers, but it is very difficult to control accurately the various parameters 
involved: positions and geometries of scatterers, values of the physical parame- 
ters; 

. numericalapproachaii3.lt allows a simpler control of the parameters and is 
fast; in practice, however, specific tools are required to perform efficient numer- 
ical computations and to render the numerical artifacts much smaller than the 
quantities of physical interest. 

The first aim of the present paper is to describe a numerical approach of this kind. 
The second aim is to highlight the efficiency of this methodology to explore the 
validity domain and limitations of analytical multiple scattering models. 

For this purpose, we will proceed as follows. In section [21 the problem of ob- 
taining random configurations is discussed; naive algorithms converge slowly when 
the scatterer concentrations are greater than 40 %. The statistical behavior of the 
configurations is determined by performing a detailed analysis of the radial dis- 
tribution function. In section [3l time-domain numerical methods are introduced. 
Elastodynamic equations are integrated using a high-order finite-difference time- 
domain scheme whose numerical artifacts are known in the case of a homogeneous 
medium. The discretization of the interfaces between the host matrix and the scat- 
terers is a key issue: special care has to be taken here to prevent the interfaces to 
introduce large numerical artifacts for physical, geometrical and numerical reasons 



111 ] . In section HI signal-processing tools are applied to the simulated data, yielding 
the effective wavenumbers. 

This numerical method is applied to a simple model of concrete, consisting of 
mortar containing composite inclusions ■ In section [H numerical experiments 
are performed at various inclusion concentrations (ranging from 3 % to 60 %), with 
longitudinal and transverse incident plane waves. Studies are performed to ensure 
that the averaged field obtained from a finite number of disordered configurations 
is representative of the theoretical effective field. In section [6l wavenumbers are 
extracted from the simulated data. Comparisons are made with multiple-scattering 
models in terms of the concentration and the adimensional frequency. In particular, 
the advantages of recent developments over the traditional Waterman- Truell model 
are confirmed. In section [71 conclusions are drawn and some future lines of research 
are described. Technical details about the computation of theoretical wavenumbers 
are given in the appendix [Xi 



2. Random configurations 
2.1. Algorithm 

Wave propagation is investigated in an infinite medium consisting of a matrix 
containing circular inclusions with a constant radius a, in the x — y plan. Matrix 
and inclusions are in perfect bonded contact; they both consist of linear elastic 
isotropic homogeneous media. 

In practice, the time-domain numerical simulations are performed in a bounded 
computational domain [Xi, X2] x [Yi, Y2]. For this purpose, N scatterers are 
introduced into the rectangular subdomain V = [Xi, X2] x [Yinf, ^sup], where 
^1 < ^nf < ^sup < ^2 (figure [T|). A minimum exclusion distance 2a + between 
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Figure 1. Computational domain. Subdomain T> containing the scatterers. A minimum exclusion distance 
2 a + g between the centers of the scatterers is assumed, where ^ depends on the mesh size. 

the centers of the scatterers is required by the numerical methods; ^ increases with 
the mesh size (section 13. ip . The x and y coordinates of the centers of circles Ci 
(z = 1, • • • , A^) are uniformly distributed in [Xi, X2] and [lint + a, Ysup — a]. Lastly, 
periodicity of the configuration is imposed along the x-axis, where the period is 
X2 — Xi. 

Two algorithms are proposed to simulate the Ci: 

Algorithm 1. 
> choose Ci randomly in T>\ 
for i = 2 to iV do 

- choose randomly Ci in T>\ 

- if CiCj >2a + ^ (j = l, — 1) then Cj is kept ; 

- otherwise choose another Cj . 

Algorithm 1 is very simple and gives quasi-uniform distributions. However, poor 
convergence is obtained at surface concentration cj) greater than 30 %, especially 
with large ^. Surface concentrations (j) greater than 50 % are beyond the reach of 



Algorithm 2. 

a compact hexagonal packing pattern consisting of N circles is 
initially introduced into T> . The side length of each hexagon is 



> repeat until sufficiently uniform distributions are obtained: 

- for z = 1 to do 

• perturb the position of Ci; 

• if CiCj > 2 a + ^ (j = 1, ■ ■ • , TV, j i^i) then d is kept . 

The mean values and variances of the Cj coordinates are measured at each itera- 
tion. When their third decimal value no long varies, then the perturbation process 
is stopped. Configurations thus-obtained close resemble to uniform distributions 
(l3l |. Algorithm 2 can be used to reach surface concentrations up to roughly 66 %. 

In practice, we recommend the use of algorithm 2 whatever the (j). In order to 
obtain a sufficiently large number of disordered patterns, the algorithm selected is 
applied M times, which gives M independent configurations for each incident wave 
at each scatterer concentration (section 12. 2p . 




2a + i; 



February 16, 2012 1:21 Waves in Random and Complex Media Wrml-Vl 



4 M. Chekroun, L. Le Marrec, B. Lombard, J. Piraux 

2.2. Radial distribution function 

In this section, we examine numerically whether the above algorithms give uni- 
form distributions. For this purpose, let us take the normalized radial distribution 
function (RDF) 

g(r) =p{r) — , 
no 

where p is the conditional probability (appendix |A]), and no is the number of 
scatterers per unit area. In an infinite statistically homogeneous domain, g 1 
as r ^ oo. The RDF is calculated numerically by counting the number n(r) of 
inclusion centers present in a circular ring with radius r and thickness Ar: 

9ir) = ^- 2a + e<r<w, (1) 

no 2 vr r Ar 

where r^ax, the distance to the nearest boundary of P, is used in order to prevent 
bounding effects. In practice, we take Ar = a / 20. The RDF is calculated for 
each inclusion in the simulation domain in order to obtain a representative value. 
In the case of dilute media, the number N of inclusions is too low to obtain a 
smooth curve, and it is not possible to determine the typical behavior of the RDF. 
The number of configurations Af is increased until a standard deviation on g{r) of 
around 5% is obtained when 7 < r/a < 10: in this range, the RDF is stabilized at 
1. The parameters of these calculations are given in table [TJ 
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Tabic 1. Parameters for the RDF calculations iQ: surface concentration (p, number of configurations J\f, number 
of inclusions N. 



As can be seen in figure [21 the RDF depends greatly on the concentration of 
the inclusions. At low concentration {(j) < 10%), the RDF can be satisfactorily 
approximated by a Heaviside function if N is sufficiently large. In other words, the 
conditional probability is uniform in this case. At higher concentrations, the local 
density of the neighbors in the vicinity of a given inclusion is increased. If i;^ < 30 %, 
this local increase is proportional to (p, as predicted by the virial expansion [l3 |. 
whereas the conditional probability is uniform if r ^ 4 a. With more densely packed 
media {(p ^ 40%), attenuated oscillations occur periodically. As the concentration 
increases, the oscillations occur farther away from the inclusion, and the period 
decreases. 

The random arrangement of circular inclusions has been studied in detail in 
many contexts. The behavior obtained here is in agreement with the results of 
statistical analyses. For further information about these distributions, see (l5| . 
The most compact 2D-arrangement corresponds to a hexagonal lattice. The RDF 
of this crystal is composed of Dirac distributions. As shown in figure [2] (right), the 
positions of the peaks in the simulated RDF do not correspond to the position of 
the Dirac distributions in the lattice. Even in the case of densely packed media, the 
heterogeneous structure cannot be approximated by a pseudo-periodic medium. 
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Figure 2. Radial distribution function g{r) at various concentrations 0, calculated with Ar = a/20 in 
JTJ. Left: RDF at </) = 6%, with Af = 3 and A/" = 100. Right: RDF at various surface concentrations <f), 
calculated with the parameters given in table [T] the thin vertical lines indicate the positions of the Dirac 
distributions in the case of a hexagonal lattice. 

3. Time-domain simulations 



3.1. Integration of elastodynamic equations 

A velocity-stress formulation of elastodynamics is followed. The physical parame- 
ters are the density p, the speeds of longitudinal waves cl and of transverse waves 



ct- The unknown are the horizontal and vertical velocity {v. 
dent components of the stress tensor (ci; 
linear hyperbolic system 



and the indepen- 
XX, (^xy, (^yy)- One has to solve the first-order 



d d d 

TT^U + A — U + B — U = 0, 

ot ox ay 



(2) 



where U = [vx, Vy, axx, o'xy, (^yy)'^ , and A and B are 5x5 matrices depending on 
the physical parameters. The system ([2]) is solved on a uniform Cartesian grid of 
Nx X Ny nodes, with mesh sizes Ax = {X2 — Xi) / Nx and Ay = (I2 — Yi)/Ny, 
and a time step At. In practice, Ax = Ay. An explicit two-step finite-difference 
ADER (Arbitrary DERivatives) scheme is used, giving fourth-order accuracy in 
both space and time (,161.1171]. With this scheme, the minimal extra distance between 
two scatterers is ^ = 3 Ax (section I2.ip . The CFL limit of stability is 



At 
Ax 



< 1, 



(3) 



where Cmax is the maximum speed of the waves in the domain. A plane wave 
analysis of this scheme has been performed in the case of a homogeneous medium 
[Ht], in terms of 9 and G = Ax/A, G €]0, 0.5], where A is the wavelength. The 
maximum artifacts are obtained when the direction of the propagation coincides 
with the grid axes, that is in the case of 1-D configurations. In this case, one has 



qie, G) = 1 



15 



1) 



4) G^ + 0(G6), 



a{9, G) 



4 7r^ 



(4) 



(02_i)(02_4^ G^ + 0(G®), 



where q is the ratio between the exact and discrete phase velocities, and a is the 
discrete attenuation [l3|- The relations (jl]) have crucial effects on the accuracy of 
the simulations, because they bound the numerical artifacts in homogeneous media. 
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3.2. Discretization of interfaces 

Three classes of drawbacks are classically associated with interfaces in finite- 
difference schemes. First, since the geometrical description of arbitrarily-shaped 
interfaces is poor, spurious diffractions are generated. Secondly, since the jump 
conditions are not enforced numerically, convergence may occur towards a non- 
physical solution. Lastly, non-smoothness of the solution across interfaces decreases 
the accuracy of the scheme, leading to spurious oscillations and even to instabili- 
ties. These three drawbacks increase with the scatterer concentration and preclude 
the use of simulations as metrological tools in highly heterogeneous media. An al- 
ternative strategy consists in using numerical methods with unstructured meshes, 
such as finite-element methods, Galerkin discontinuous methods, and spectral ele- 
ment methods @, [l^. But the computational cost of these methods would be much 
higher due to the meshing, and the stability condition ([3]) is penalized. 
To overcome these drawbacks, we use a r-th order immersed interface method 



ll|, |20(]. This numerical method modifies the ADER scheme at grid points close to 
the interfaces, based on the jump conditions up to the r-th order, the elastodynamic 
equations, and the Beltrami equations. This procedure associates the efficiency of 
Cartesian grid methods and the accuracy of an interface meshing. The work is 
mainly carried out during a preprocessing step, before the numerical integration. 
At each time step, 0{C / Ax) matrix-vector products are done, where C is the 
total perimeter of the interfaces, and the matrices are small, typically 5 x 100. The 
results are then injected into the scheme. After optimizing the codes, the additional 
CPU time required by the immersed interface method can be made negligible in 
comparison with the CPU time required by the scheme (less than 1%). 



3.3. Intensive computing 

To obtain reliable effective wavenumbers values, the numerical methods used must 
meet the following specifications: 

• large computational domains, such as grids consisting of 10^ x 10^ nodes, involv- 
ing 10 Go of data; 

• long integration times, consisting for example of 10^ time steps; 

• several simulations, so as to increase the number of independent disordered con- 
figurations, to TV = 3, for example; 

• processing a large number of scatterers with the immersed interface method: 
1500 interfaces when (j) = 48 %, for example; 

• performing many simulations in the parametric studies, e.g. in terms of the 
scatterer concentration or the incident wave polarizations. 

To meet these specifications, the computer codes are parallelized. Domain de- 
composition is performed in the x-direction, associating each slide with one com- 
putational process. All the slices have the same size and contain approximately 
the same number of inclusions, so that the computational cost of each process is 
roughly the same. After each time step, data are exchanged between neighboring 
processes. 

In practice, the simulations presented in sections [5] and [6] were performed on a 
cluster of 4 PC bi-processor quadricores, amounting to 32 processes. The optimum 
speed-up 32 was reached. The communication time between processes was negligi- 
ble in comparison with the computational cost of each process. After parallelizing, 
the configurations investigated in section [5] required 1 hour of preprocessing (due 
to the use of the immersed interface method) and 24 hours of integration (due to 
the use of the ADER scheme). 



February 16, 2012 1:21 Waves in Random and Complex Media Wrml-Vl 

Waves in Random and Complex Media 7 

4. Data processing 

4.1. Numerical coherent field 

At each time step, the components of U have to be stored inside the subdomain 
containing the inclusions. For this purpose, a uniform network consisting of A'^; 
hues and Nc columns of receivers is placed in the subdomain D. The position of 
the receivers is given by (x, = Xi + i Ac, Uj = Ylnf + j A;), where i = 0, . . . , A'c — 1 
and j = . . . , Ni — \. The bottom and top receivers are sufficiently far from the 
sides of the computational domain to prevent spurious effects from being recorded. 
These columns of receivers are visible in figure [5]- (a). 

The acquisition setup has to meet some specifications in order to prevent the 
occurrence of aliasing and low resolution problems [21]. Aliasing occurs when the 
distance A/ is larger than the shortest wavelength under consideration, while the 
resolution is limited by the total length Ni A; of the acquisition setup (i.e. the 
distance between the first and last receivers). 

The field recorded on each array (each column of receivers) corresponds to a 
field propagating in a given disordered configuration. Summing the time histories 
of these Nc arrays in all the M simulations gives a coherent field propagating in 
the y direction. The relationship between the coherent field and the effective field 
will be discussed in section 15. 3i 

The polarization of the coherent field is the same as the polarization of the 
incident plane wave. In the following section, we will therefore deal with either the 
L (longitudinal) or T (transverse) case alone. 

4.2. Extraction of the coherent wavenumbers 

The coherent phase velocity c{u!) is computed by applying a p — a; transform to 
the space-time data on the coherent field, where p is the slowness of the waves 
(p = 1 / c) and CO is the angular frequency [12, H^] . The time Fourier transform of 
the coherent field s{yj, u) is denoted by 

s{y„co) = ^(y„.;)e-*-P''(-)^% (5) 

where A{yj, co) is the amplitude spectrum at yj, and po(w) needs to be determined. 
A p— a; stack quantity s(p, oj) is then defined as 

Ni-l 

s{p,u;)= J]^(y„t^)e*-(P-P°(-))^^ (6) 

j=0 

s(p, (jj) is computed at several p values. Given w, the maximum value of the modulus 
|s(p, a;)| is reached at p = po('^) = 1 / c(a;). The phase velocity dispersion curve is 
then obtained by taking the maximum locus on the 2-D map |,s(p, a;)|. An error 
estimate is also deduced [13]. 

Let us now examine the attenuation of the coherent field. This quantity is esti- 
mated from the decrease in the amplitude spectrum of the coherent field during 
the propagation of the waves. In the frequency domain, the amplitude A{yj, u) in 
©-([G]) is assumed to satisfy an exponential decay with distance 

^(y„u;) = Ao(o;)e-"('^)(^^-^"), (7) 



where Aq(uj) is the amplitude of s at the first receiver located at the offset yo, and 
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a{uj) is the attenuation factor. Although two points suffice to be able to calculate 
Q, the slope of a least-squares fit of ln(A(yj, u)) over the whole range of reception 
gives more accurate results. An error estimate is also deduced. 



5. Numerical experiments 
5.1. Validation 




Figure 3. Validation test 1/2. Numerical dispersion (a) and numerical attenuation (b) in a ID homoge- 
neous medium: analytical values l(4]| in red lines, and measured values in blue circles. 



In the first test, transverse wave propagation was simulated in a homogeneous 
1-D cement matrix. The dispersion and attenuation measured were due only to 
numerical artifacts occurring in the ADER scheme. Comparisons between the the- 
oretical (jlj) and measured dispersion and attenuation values is made in figure El 
The error between the theoretical and measured curves is less than 10~^% in the 
frequency range of interest. The signal processing tools used and the acquisition 
setup chosen are therefore suitable for accurately assessing the dispersion and the 
attenuation, and the risk of adding significant signal processing artifacts is thus 
avoided. 




-0.04 -0.02 0.02 0.04 1.2E-5 1.6E-5 2E-5 2.4E-5 2.8E-5 



x(m) t(s) 



Figure 4. Validation test 2/2. Snapshot and time-history of axx'- comparison between the numerical and 
the analytical solutions at the receiver R. In (a), the green-red palette and yellow-magenta palette denote 
L waves and T waves, respectively. 



In the second test, the wave propagation was simulated in a 2-D medium with 
a single scatterer centered at (0, 0). The source was the plane compressional 
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wave described in section [5l The diffracted fields were stored in tlie receiver R 
at (—0.25, —0.25). Figure H] gives a snapsliot of tlie stress cr^x after 600 time steps 

(a) and compares tlie numerical and analytical values of axx during 1400 time steps 

(b) . The exact solution was computed by performing standard Fourier-Bessel de- 
compositions. The excellent agreement observed confirms the validity of both the 
ADER scheme and the immersed interface method. 



5.2. Numerical setup 

The numerical method presented in sections [2] to S] is now applied to some physi- 
cally relevant configurations. A simple model of concrete is studied, where circular 
aggregates with a radius a = 6 mm are embedded in a homogeneous cement matrix. 
The physical parameters are 

{po, co,Lj co,t) = (2050 kg. m^'^, 3950 m.s^"*^, 2250 m.s^"*^) in the cement matrix, 
(pi) ci,L, ci,t) = (2610kg. m~^, 4300 m.s^"*^, 2470 m.s^"*^) in the aggregates. 

A parametric study is performed in terms of the concentration, from (p = 3% to 
60%. The domain of investigation presented in table [2] is discretized on Nx x Ny = 
7200 X 7200 nodes, hence Ax = Ay = W'^ m. The CFL number ([3]) is 6* = 0.95, 
giving At = 2.21 10~^ s. A third-order immersed interface method is implemented 
(r = 3 in section [3^ . The source is a plane longitudinal (L) or transverse (T) wave 
propagating along the y-axis, initially outside the domain T> (figures [1] and [5]-a) . 
The time evolution of the source is a Ricker with a central frequency 250 kHz. The 
frequencies of interest range between 50 kHz and 600 kHz. 



-0.36 -hO.36 -0.02 -^0.7 +0.2 +0.6 

Tabic 2. Coordinates of the physieal domain and those of the subdomain X> (scction[2j, in meters. 



The acquisition network contains = 100 columns and Ni = 400 lines, with 
the spacing = 0.001 m and Ac = 0.072 m, respectively. Each column is an 
array of receivers that follows the wave propagation in a particular realization of 
disorder. Performing M = 2> simulations yields 3 x 100 independent disordered 
configurations. This acquisition setup gives the following bounds on the standard 
errors: from 2m.s~^ at / = 50 kHz to 0.2m.s~^ at / = 600 kHz in the case of the 
phase velocity, and from 0.05Np.m~^ at 50 kHz to O.lNp.m"^ at 600 kHz in that 
of the attenuation. 

Signals recorded along 2 different arrays logically show different behaviors. Figure 
[6]-a gives the time histories at various receivers along one particular array. A main 
wave train is clearly visible in each of the time histories, followed by a coda. 

As recalled in section 14.11 a coherent signal can be obtained by averaging the 
signals recorded on the various arrays 25||. The Vy component is used in the case 
of an incident L-wave, whereas the Vx component is used in that of an incident 
T-wave (stress components provide the same results). An example of coherent 
seismogram is presented in figure [6]-b. The coda has disappeared, and the main 
wave train behaves like a plane wave propagating in a homogeneous (but dispersive 
and attenuating) medium. 
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(b) 





Figure 5. Incident L-wave, concentration </> = 42 %, initial instant (a) and after 3000 time steps (b). The 
vertical columns denote the positions of the receivers. 
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Figure 6. Incident plane T-wave in a medium with 24 % inclusion concentration: (a) signals recorded along 
an array, (b) coherent signals obtained after summation. 

5.3. Convergence of the coherent field to the effective field 

The coherent field is obtained by averaging the signals recorded along the 300 
arrays of receivers. If this number is too low, the estimated properties will still 
be dependent on the configuration of the scatterers encountered. Theoretically, the 
effective wavenumber can be defined by taking an infinite number of configurations, 
which means that all the possible configurations of scatterers will be taken into 
account; but this approach is obviously impracticable. The aim of this paragraph 
is to show that 300 arrays suffice for estimating the effective field and hence, the 
effective wavenumber. 

In each case (in terms of the density (p of the scatterer and the incident wave 
/3 = L,T), a coherent signal is computed with an increasing number Na of arrays 
ranging from 1 to 300. The Na arrays are chosen randomly among the 300 available 
ones, to avoid taking consecutive arrays which are located too near each other in 
the medium. The properties and a/j are then evaluated from the averaged signal, 
and their evolution with Na is then studied at a given frequency. Figure [7] shows 
how aT evolves with Na, at = 36% and / = 300 kHz. As this evolution depends 
on the Na arrays selected in the averaging procedure, the study is repeated 10 
times. 





+0.2 Np.m ' ' 




a -0.2 Np.m"' 

e 



50 100 150 200 250 300 

N 



Figure 7. Evolution of ot in terms of Af^ at / = 300 kHz, tj) = 36%. 



The value obtained by summing the fields over the 300 arrays is taken as a 
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reference value for the effective medium: in this case, it was ae,r=3.06 Np.m^^. In 
figure [71 the red curve gives the envelope of all 10 curves in gray, corresponding to 
2 standard deviations of aT^Na). As was to be expected, as Na increases, the value 
of ar tends towards the reference value ae,T- This figure shows that when there 
are too few configurations, ae,T cannot be accurately assessed; for instance, taking 
only 50 configurations results in an uncertainty greater than INp.m"^. Taking 
±0.2 Np.m^^ to be an acceptable level of uncertainty for Oe^T, the optimum number 
Nopt of configurations requested must be greater than Nopt — 210 in the case of 
the present example. 
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Figure 8. Optimum number Nopt of eonfiguration of seatterers for determining or at all densities. 



Figure [8] summarizes the values of Nopt obtained in all the cases studied, at 
several frequencies, using the procedure above. The criterion used to obtain an 
accuracy of about Ce,/3 was ztSm.s""*^; in the case of ae,/3) we took ±0.2 Np.m""^. 
Except for the higher frequencies (> 400 kHz), only about 20 configurations are 
required for assessing the effective velocity phase Ce,/3, whereas greater values of 
Nopt are required for assessing the effective attenuation. If Nopt is obtained at a 
given scatterer concentration, its value will increase with the frequency at a given 
polarity, and it will be almost twice as high with T-waves as with L- waves at a given 
frequency. Convergence therefore depends mainly on the size of the wavelength, as 
the effects of multiple scattering are greater at shorter wavelengths. 

The results obtained on ot at high frequencies (500 and 600 kHz) show that 
the optimum number of configurations Nopt was almost 300, which was the max- 
imum number of configurations available, whatever the density of the seatterers. 
In the present 2 cases, it was not possible to say whether the optimum number 
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of configurations was actually reached, so as to be able to assess the attenuation 
with a sufficiently high level of certainty. In all the other cases, the Ce,/3 and ae,/3 
values obtained based on 300 scatterer configurations were fully representative of 
the effective wavenumber. 



6. Numerical results 

6.1. Simulated effective wavenumbers 

Results obtained in the numerical simulations with the two polarizations of the 
incident wave (/3 = L, T) and at various scatterer concentrations (f) are presented 
in figure [H In a first approximation, the effective phase velocity was found to 
be proportional to the density of the inclusions, increasing monotonically with 
(p. The effective phase velocity also showed a dispersive behavior, which became 
more conspicuous as (p increased. This effect was stronger at the lower frequencies 
(at fco,La ^ 1 and /co,t« ^ 2), where the value of the overall phase velocity was 
lower than at high frequencies. A maximum value of Ce was reached at ko^ia ~ 1 
and kQ^TCL ~ 2. At higher frequencies, Ce,/3 remained almost constant but small 
fiuctuations are visible, up to lOm.s"^ at 60% which amount to less than 0.5% 
of the phase velocity. However, the positions of these local extrema are almost 4>- 
invariant. At high frequencies, the mean phase velocity c = (1 — ^) cq,/? + 4>ci^i3 is a 
good approximation in the case of dilute media and an upper limit in that of more 
densely packed media. 

The attenuation is more difficult to explain: contrary to what occurs with the 
phase velocity, the attenuation does not depend monotonically on the concentra- 
tion. The frequency dependence shows 0-invariant peaks corresponding to the local 
maxima of the phase velocity, mainly at ko^ia ~ 1 and kQxa ~ 3 in the case of L- 
waves, and /cq^to ~ 2 and /co,Ta ~ 5 in that of T waves. In dilute media (cp < 20%), 
the attenuation is proportional to (p. In denser heterogeneous media, the behavior 
depends on both the inclusion concentration and the frequency range. Three types 
of overall behaviors were observed: 

• around the previously mentioned peaks, ae,/3 increases with (p; 

• at high frequencies between these peaks, ae,/3 remains at an almost constant 
value when 30% < (p < 60%; 

• at low frequencies, ae,^ reaches a maximum at w 30%, and then it decreases. 

With an incident T-wave and fco,r a > 7, the behavior is less clear-cut: the attenu- 
ation reaches a peak at around (p — 30%. However, the accuracy of these findings 
was not confirmed in the section 15.31 

6.2. Comparison between various theoretical models 

In this section, the effective wavenumbers predicted by the Independent Scattering 
Approximation (ISA), Waterman- Truell (WT) [3], and recent Linton-Martin (LM) 
[3] multiple-scattering models are presented. Technical details can be found in the 
appendix El The results obtained with these models are compared in figure [TOl in 
the case of a longitudinal incident wave, and the same comments apply in the case 
of transverse waves. 

ISA and WT models give similar results. The phase velocities differ only at low 
frequencies in the case of densely packed media, where both models are inaccurate; 
we will therefore focus on the WT predictions. As regards the phase velocity, all 
these models predicted the same overall behavior as the simulations: a monotonic 
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Figure 9. Effective properties obtained with numerical simulations and signal processing tools, with lon- 
gitudinal (left) and transverse (right) incident waves at various inclusion concentrations </>. Top: phase 
velocity the horizontal dashed lines give the mean phase velocity c. Bottom: attenuation ae,/3 where 

</> = 30% (that is, around the critical threshold mentioned in the discussion) is given by a dashed line; the 
insert is a zoom in the low-frequency range. 



increase with the concentration, the same i;^-invariant position of the local extrema, 
and an increase in the dispersion with the concentration. However, the attenuation 
given by the ISA and WT models are mainly linear functions of the concentra- 
tion, contrary to what observed with the simulations: the attenuation occurring 
in densely packed media was clearly over-estimated. The attenuation predicted by 
the LM model differed considerably from that obtained with the previous mod- 
els. First, the attenuation showed a linear dependence on the concentration only 
at low densities ^ < 30%, where it reached a maximum (except for the peak at 
^o,Lfl — !)• Secondly, the attenuation was negative in some (low or high) frequency 
ranges, which is unphysical, but this occurs only at very high concentrations. 



6.3. Comparison between theoretical models and numerical results 

These comparisons between models and simulations made it possible to determine 
the range of validity of each model (see figure [TT] for the phase velocity and figure 
[12] for the attenuation). In view of the signal processing limitations, satisfactorily 
results were defined as those with an error of less than 5m.s~^ in the case of the 
phase velocity, and less than 0.2Np.m~^ in that of the attenuation. On this basis, 
the WT model is suitable for dealing only with very dilute media (j) ^ 12% for Cg, 
and (p ^ 6% for Og, whereas the LM model gives accurate results up to inclusion 
concentrations of 24% in terms of both the phase velocity and the attenuation. This 
considerable difference is probably attributable to the hole correction occurring in 
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Figure 10. Effective properties of the longitudinal ineident waves obtained writli the Independant Scatter- 
ing Approximation (ISA), Waterman- Truell (WT) and Linton-Martin (LM) models. 
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Figure 11. Comparison between numerical simulations and modeling predictions of the phase velocity in 
the case of longitudinal (left) and transverse (right) incident waves at the lower inclusion concentrations 
(</> < 30%). 
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Figure 12. Comparison between the effective attenuation obtained in numerical simulation and modeling 
predictions at the lower inclusion concentrations {(p < 30%). Results obtained in the case of longitudinal 
(left) and transverse (right) incident waves versus the Waterman- Truell (top) and Linton-Martin predic- 
tions (bottom). Solid line: results of numerical simulations; dashed line: modeling predictions. 

the LM model. As mentioned above, the simulation at </> = 6% was performed using 
only three realizations of the simulation domain, which resulted in the non-smooth 
RDF shown in figure [21 however, the LM and WT models both gave excellent 
results at = 6%. This indicates that the hole correction has the most significant 
effects Sit (j) > 10%. At (p > 24%, there is a marked discrepancy between the RDF 
and the Heaviside step function. A more realistic form of p(r2|ri) presented in ()Aip 
might extend the range of validity of the LM model to include higher densities. 
The WT and LM models are often referred to second-order models, because of the 
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Figure 13. Error (/3 = I/,T) versus iji at various adimensional frequencies (top: WT; bottom: LM; left: 
longitudinal wave; right: transverse wave). X: implicit method; o: explicit method. 



second-order Taylor expansion (1A5[) used to estimate k^.^^. However, to our knowl- 
edge, this definition has never been justified either theoretically or experimentally. 
We will therefore examine this question numerically. Let us assume 



(8) 



where fce,^ is the effective wavenumber obtained in the numerical simulations. The 
difference between the numerical and the theoretical results are defined by 



||(fce,/3a) - (fee,/? a) II, 

o^lho (5i,/3 - C^l,/3) + (?^o)^ ('^2,/3 - ^2,/?) || + 0(n\), 



(9) 



where fce,/3 is the wavenumber deduced from the LM or WT model. This difference 
is plotted at various (adimensional) frequencies versus the surface ratio cp = uq -no? 
in figure [131 With the WT model, Ep is governed throughout the whole frequency 
and concentration range by a slope of 2 decades per decade on a log-log scale: 
therefore depends mainly on 52,/3 — ^^2,/3- In other words, the WT model predicts d2,/3 
inaccurately, even at low densities. In fact, there is no reason why the WT model 
should be preferred to the ISA model, which is only a first-order model. With the 
LM model, £p is independent of the inclusion concentration up to i?!) ~ 18%, and 
has a -|-3 slope on log-log scale. This confirms the accuracy of the second-order 
coefficient 1^2, /3 obtained with the hole correction. 
When Efi was calculated using an implicit formulation to obtain the wavenumber 
[3]) . was found to have the same concentration-dependent properties with both 
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the WT and the LM models (x -curves in figure [T3]l . The imphcit formulation did 
not improve the quality of the results, and LM is therefore an intrinsically-second- 
order model. 



7. Conclusion 

The effective properties of random elastic media were calculated here using purely 
numerical methods. Combining sophisticated methods of simulation (the fourth- 
order ADER scheme and the immersed interface method) and signal processing 
tools yielded reference solutions for both the real and imaginary parts of the ef- 
fective wavenumbers. With this approach, the accuracy of the simulations does 
not depend on the scatterer concentration. Maximum computational efficiency is 
obtained by performing domain decomposition and parallelizing the algorithms. 

This numerical method was applied in the present study to a 2-D model of con- 
crete. In this case, the numerical simulations confirmed that traditional models 
(such as the Waterman- Truell model) are valid roughly up to inclusion densities 
of 10%, whereas the recent Linton-Martin model (extended by Conoir-Norris to 
elastodynamics) is valid up to densities of 25 %. In particular, the present simula- 
tions confirmed the validity of the second-order term in the Linton-Martin model, 
as previously done in a theoretical study 0] ■ 

The numerical method presented here can be used to handle more complex con- 
figurations, such as a granulometry or composite containing scatterers of various 
shapes. Other constitutive laws could also be introduced, such as viscoelastic laws 
accounting for dissipative effects [1^. Lastly, the possibility of extending the present 
approach to 3-D configurations is a great computational challenge. Preliminary 
tests have already been conducted on the numerical methods with fiuid scatterers 
included in a fiuid matrix. 
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Appendix A. Effective wavenumbers 

As established in the scattered field ips recorded at r and due to an inclusion 

centered at ri can be expressed through a linear operator T applied to the excit- 
ing field ijje'- ^s(r;ri) = ^(ri) i/;e(r;ri). If the positions of the N inclusions are 
known, closed-form solution of the problem can be obtained. If, on the contrary, 
the inclusions are randomly distributed, it is generally attempted to determine the 
effective field (V'e) corresponding to the ensemble average of the positions of all the 
inclusions. The effective field at one representative inclusion (say the first one) is 
expressed in terms of the scattering induced by another representative inclusion 
(say the second). Based on the quasicristalline approximation (QCA) the lat- 
ter scattered field is assumed to be excited by the same effective field as the first 
inclusion: 

(V'e(r;ri)) =V'*(r) + (iV-l) / r(r2) (V'e(r; rs)) p(r2|ri) dra, (Al) 
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where "01 is the incident field. The probabihty density p(r2|ri) in (lAlh . which is 
a pair-correlation function, expresses the probability of finding an inclusion at r2, 
given that an inclusion is placed at ri. 

Various models have been developed for determining the effective wavenumbers 
analytically, such as the Waterman- Truell (WT) pj] and Lloyd-Berry (LB) (29j | 
models, to cite but a few. As pointed out by Linton-Martin (LM) [3], these mod- 
els differ mainly in the hole correction, i.e. the assumption made about p in the 
integration of (lAip . The LB model is based on the following assumption: 

if |r2 - ril <b, 

b = 2a + ^, (A2) 

1 if |r2 — ri| > 6, 

where a and ^ are defined in section 12.11 To obtain the effective wavenumbers in L 
and T waves, the effective field can be decomposed into a modal sum, where the 
modal amplitude depends on the effective wavenumber 0]. Introducing this form 
into ()Aip and using (IA2p . the problem reduces to searching for the non-trivial 
solution of the infinite linear system 0, equation (31)]: 



iV- 1 
no 



P(r2|ri) 



det (I - 2 no M T) = 0, (A3) 
where I is the identity matrix, T is a matrix defined in P, equation (37)], and 



M 



Ml 
Mr 



[m,n] — ^2 _ y^eb J'„i_n(.^eb) H'^^_^{kfsb) — kj^b Jm-n{keb) hU_^' {kpb)j . 

(A4) 

In practice, the modal sum in (|A3p is truncated. The effective wavenumbers ke 
satisfying ()A3|1 are associated with ke^L and ke^x- However, searching /ce,L and ke,T 
is an intricate and time-consuming process. Another explicit but approximate form 
can be obtained using Taylor expansions: 

^e,/3 = ^1 + "-O^l,/? + ("-o)^ '52,/3, (A5) 

where 5i^p and 82,13 are defined in [Q, equations (62a) and (62b)] in the LB model. 
Let us now consider the WT model, which is based on 

iV - 1 , , , [0 if (r2 -ri) - ey < r/, 

p(r2|ri) = < with r/ — >■ 0. (A6) 

^0 [ 1 if (r2 - ri) - By > r], 



where (r2 — ri) • ey is the distance between two inclusions in the direction of 
propagation (ey here) of the incident wave. Based on this hypothesis, the system 
will have the same form as ([A3]), but M will be different, and this parameter is 
deduced from equations (12-13)]: 

M ^ ( ^ (-1)'"""\ ^ 

^ ' ikg \ke - kg ke + kg J ' 



The approximate form (jA5|) can still be used with the WT model. The same ex- 
pression for 61, g in (jASP can be used as with the LB model, but the definition for 
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(52,/3 in this case is that given in [4,, equation (4)]. Note that ISA matches with (|A5p 
to the first order: fe^ ^ = fe^ + no 5i^p. 
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